All the adiabatic bound states of NO2 (J=0) 
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£N| ■ We calculated all 2967 even and odd bound states of the adiabatic ground 

o, 

state of NO2 , using a modification of the ah initio potential energy surface 
q^ ■ of Leonardi et al. [J. Chem. Phys. 105, 9051 (1996)]. The calculation 



cr 



X 



was performed by harmonic inversion of the Chebyshev correlation function 



generated by a DVR Hamiltonian in Radau coordinates. The relative error 
for the computed eigenenergies (measured from the potential minimum), is 



10 or better, corresponding to an absolute error of less than about 2.5 
cm -1 . Near the dissociation threshold the average density of states is about 
0.2/cm _1 . Statistical analysis of the states shows some interesting structure 
of the rigidity parameter A3 as a function of energy. 
Pacs 34.10.+X, 34.30. +h, 05.45.+b 



* Author to whom correspondence should be addressed 



I. INTRODUCTION 

The complexity of the spectrum of NO2 has always been a challenge, both to experimen- 
talists and theorists. Well-resolved spectra were first measured only after the advent of the 
laser, and the introduction of cooled molecular beams allowed the measurement of complete 
spectra in certain energy ranges (still not up to the dissociation limit) JI|-Q. One reason 
for the spectral congestion is the conical intersection (found by Gillespie et al. in 1975 [BJ) 
mixing the X 2 A\ electronic ground state surface at an energy of only 1.289 eV with the first 
excited A 2 B 2 surface. The dissociation energy Vdi SS of X 2 A\ is 3.226 eV. In the experiment, 
this leads to disturbed progressions and increased state density, while the theorist is forced 
to include the vibronic coupling in the calculation of more than the lowest few vibrational 
states of the molecule |||7|]. 

There are only a few molecules for which anything like a full list of vibrational states (even for 
J = 0) has been calculated. Examples are the pioneering work of Tennyson et al. || on H3", 
and the more recent calculations on HO2 0, HCN |TO| , H^ JT1J and NO2 fL2 |. The reason 



for this situation is, of course, the immense computational expense of such calculations, 
which one generally applies only to potential energy surfaces (PES) which are quite well 
known. These, however, are still scarce. 

In this paper we present first results of an effort to compute bound and resonance states 
of NO2 • Two developments were of great help to us: first the calculation of good PESs 
for the two lowest electronic states of NO2 including their coupling by the Siena-Wuppertal 
group fl2l1 ; second, the advent of the filter diagonalization method [|13HT7H for solving large 
eigenvalue problems. The first step of this endavour, which we report here, is the calculation 
of all (~3000) adiabatic bound states of the ground PES. This in itself is non-trivial, since 
near the dissociation limit the average state density is ~0.2 states per cm -1 for each of the two 
symmetry species. The computation of non-adiabatic energy eigenvalues, which includes the 



surface coupling, is under way JT8[ , and the method itself allows the introduction of imaginary 



absorption potentials, and, therefore, the extension of the calculations to resonances. (This 



has been tested for H3 [jiTf ) 

Another benefit of this project is the possibility to perform a thorough statistical analysis 

over the whole energy range of a molecular system. This has so far been almost impossible 

for real molecules, because of missing data in the experiments, or an insufficient number of 

eigenstates in the calculation, both leading to large statistical errors 0,[|||[7|. 

In the following we discuss first some properties and necessary modifications to the PES of 



ref. ||12| , then the computation, and, finally, some results on the statistics of the computed 



levels. 



II. THE POTENTIAL ENERGY SURFACES 

Since the two lowest adiabatic PESs of NO2 have a conical intersection at an energy of 
about one third of the dissociation energy, they must both be considered in any calculation 
of more than the lowest few vibrational states |J. However, the intersection makes their 
direct interpolation impractical. Instead, one interpolates the diabatic surfaces Vn and V22, 
and, in addition, their coupling V12. From such fits one computes the energies of the two 
adiabatic surfaces Vx (ground state) and V^ (first excited state) whenever they are needed 
from 



V^ff = \{V n + V 22 ) T f-(Vn - V 22 y + V? 2 (1) 

The surfaces used in our calculations have been derived in several steps from nearly 1000 ab 
initio points computed in ref. |L9||. These points have been diabatized [|2Q| , [T~9|1 by a unitary 



transform, demanding that certain dipole matrix elements be small. All computed points are 
in the region (here we use bond coordinates) 2.08ao < ri, r-i < 2.50ao, 70° < (3 < 180°, called 
region D in ref. [^OJ. To produce surfaces useful for the calculation of vibrational levels, they 
had to be fitted to an analytical form, furthermore they had to be extrapolated to regions 
outside of region D, i.e. to shorter and larger interatomic distances, and to small angles. 
This has been done in several iterations by the Siena and Wuppertal groups [|I^,[n],[EJ, who 



also utilized empirical data from the measured NO2 spectrum (see ref. [1] in |TI|) to improve 
the final fits of ref. JT2] . It is these fits, called here V^ fc p (for Leonardi-Petrongolo) from which 
we started. 



While the authors of ref. [jE^ had no problems to compute vibrational energies with their 
fits, early test calculations with our method showed that the short- and long-range extrapo- 
lations in V^f have (unphysical) kinks, which cannot be tolerated when a discrete variable 
representation (DVR) is used. Therefore, we modified the surfaces by smoothing out the 



kinks as described below. Preliminary calculations done by C. Petrongolo |22|] with this 



surface using the code from ref. |12|| , and our nonadiabatic computation |L8| show an equally 



good agreement with experimental data at energies up to 1.22 eV. Additionally, our surface 
modification improves the computed state density up to 2.3 eV, which had been too high 
compared with the experimental one in ref. ||12|| . 

In detail, the following modifications were applied (energies are hartree units = 27.211396 
eV, distances bohr units a = 0.529177-10~ 10 m): 

(a) The short-range correction (for r, < 2.08) of the lower diabatic surface V{f was changed 
to (shown only for r 1; but similar for r 2 ): 

V 11 (r 1 ,r 2 ,P) = V£ p (r 1 ,r 2 ,/3){l-Q(a,bri)) (2) 

+ (6.00(n - 2.08) 2 - Fn P (2.08, r 2 , (3))Q(a, bn) 
V£ P E D, 

here the function Q(a, br) = 1 — P(a, br) is the incomplete gamma function cf. ref. p3| ) and 
the parameters a and b are listed in Table [p. 

(b) The short-range correction to V22 was dropped, and V 22 (ri,r 2 , ft) set to a finite high 
value V 22 (l.50,r 2 , j3) for r\ < 1.50, and similarly for r 2 < 1.50. 

(c) The small angle (/3 < 70°) corrections to the potentials Vn and V 2 2 were dropped 
completely. 

(d) For small and large angles (/3 < 70° and (3 > 131°) the coupling potential was smoothly 
set to zero: 



Vni/3) = ^| p (/3)(l - Q(a, 6(/3 - 71°)))for (3 < 71° (3) 

V 12 {{3) = ^ p (/3)(l - Q(a, 6(/3 - 131°)))for (3 > 131°. (4) 

The values used for the parameters a and b are listed in Table |. 

(e) The long-range correction of the diagonal potentials V„, which are necessary to 
produce the correct dissociation limits, where modified by changing the switching function, 
now reading (for n > r max = 3.00) 

V tl (r u r 2 ) = (1 - Q(l, b(n - r max ))^P + 

Q(l, 6(n - r max )(V diss + V N0 (r 2 )). (5) 

and similar for r 2 . The parameters a and b are again listed in table |. 

III. COMPUTATION 

For computing bound states of N0 2 we use the low-storage version of the filter- 
diagonalization method introduced recently by Mandelshtam and Taylor [16]. This method 



is conceptually based on the filter-diagonalization procedure of Wall and Neuhauser WA 



which extracts the system eigenenergies by harmonic inversion of a time correlation function 
C(t). The method of ref. [IB| is designed for a direct harmonic inversion of the Chebyshev 
correlation function |24[] 

c n = {€o\T n (H)\£ ) ~ ^2 d k cos nu k , (6) 

k 

for the eigenenergies E k = coso^ and amplitudes d k - The computation of the c n sequence 
is done to essentially the machine precision using a very inexpensive iterative numerical 
scheme, 

6 = H£ , £ n+ i = 2H£ n - £ n _i , (7) 

with c n being generated using c 2n = 2(^„|^„) - c , c 2n _i = 2(^ n _i|^ n ) - c x . This requires 
to store only a few vectors at a time, assuming that the matrix-vector multiplication is 



implemented without explicit storage of the Hamiltonian matrix. The spectral analysis part 
(i.e., the harmonic inversion of c n by the filter-diagonalization) is carried out independently 
and computationally very efficiently after the sequence c n is computed. All these features 
assure the very high performance of the overall numerical procedure. 
For the Hamiltonian matrix representation we choose Radau coordinates (Ri, R 2 , G) in which 



all mixed derivatives in the kinetic energy vanish [Pq] , to account for the C 2v symmetry of 
the NO2 molecule. The fast application of the Hamiltonian matrix to a vector is achieved 
by implementing a potential adapted DVR, which we set up in a way which allowed us 
to identify a single parameter, the potential cutoff energy V cu t, on which all other grid 



parameters depend. The bases for the radial parts of the Hamiltonian are sinc-DVRs |26 
characterized by their spatial extensions Rimin to Ri max (i = 1, 2), and sizes n\ = n 2 , which 
are diagonalized and then truncated by an energy cutoff to n lb = n 2b . For the angular part 



we use a Legendre-DVR |27] characterized by its size n^. The angular extension is always 
0°to 180°. 

Convergence studies in which all grid parameters are varied independently are prohibited 
in this study by the computing times involved. Instead, these parameters were fixed in ID 
calculations, but in a consistent way: We first set an accuracy goal, in this case a relative 
accuracy of 10~ 4 for states just below the true NO2 dissociation energy of 3.226 eV. With 
this desired accuracy we adjust (in ID calculations at the equilibrium angle Bo) the spatial 
extensions Ri m in to Rimax (i = 1, 2) of the radial grids, while the sizes rij are kept very large. 
Then the sizes of the spatial grids, ni = n 2 , are adjusted to yield this same desired accuracy, 
similarly for the angular grid (at R\ = R 2 = Ro), whose extension, 0°to 180°, is kept. The 
equidistant grid spacings of the 1D-DVR points AR, AG, on the other hand, determine a 
maximum kinetic energy, 



7V 2 
-radial 



T — 2MAW and (8) 



n 2 ( 1 1\1 



^angular ~ l | (q\ 

max ~2fi\R 2 , R\) AG 2 ' K) 

(We used R\ = Rq and R 2 = 00 as an estimate for the moment of inertia in eq. |9], which 



seems reasonable for an almost dissociating wavepacket.) 

The larger one of these kinetic energies, which in our example (10 -4 relative accuracy at 
3.226 eV) are 4.5 eV for the radial grids and 8.0 eV for the angular grid, is then used to 
fix the potential cutoff V cu t- For consistency, the radial grids are increased in size such that 
T^ 1 = T^f lar . Then all DVR-points for which V(R U R 2 , Q) > V cut are deleted. As a 
final result for the accuracy goal mentioned above we get a primitive grid of m = n 2 = 127, 
Rimin = Rimax = 1.5a , R\ ma x = Rimax = 5.5a , n 3 = 275, a ID truncation of the radial 
grids to nib = n 2 b = 100, and a final truncated grid of 1 112 707 points. 
A further parameter of the method is the number of iterations of the filter diagonalization, 
eq. [7|, which is one half of the number of Chebyshev coefficients c n . In addition, for the final 
stage of the calculation the window parameters E min , E max , and the number of window basis 
functions (described in ref. ||28| ) have to be fixed. This was done such that the extracted 
eigenenergies were converged to 6 figures with respect to the sequence of coefficients c n , i.e. 
two orders of magnitude better than our overall accuracy goal. This led to a number of 
300 000 coefficients, whose calculation for the 8 eV set of one symmetry needed about 30 
days of CPU-time on an IBM RS6000/59H. 

To ensure that our overall accuracy is in the exponential convergence region, and not in some 
- misleading — small local convergence minimum |2J| , we ran four different grids consistent 



with energy cutoffs of 4, 6, 8 and 8.1 eV. The last one has been used in comparison to the last 
but one to observe numerical noise, and to have an independent set of eigenvalues ensuring 
that we do not miss any levels. All inconsistencies between the two largest sets have been 
cleared. 

The errors of the computed eigenvalues can be approximately taken as the noisy difference 
between the calculation with V cut = 8.0 eV and the one with V cut = 8.1 eV. As Fig. [J] shows, 
it increases with energy above about 2 eV, and above 3 eV has a plateau of about 10 -4 eV. 
This plateau, which is also shown by the differences of the smaller calculations (with V cu t 
= 4.0 and 6.0 eV, respectively), is caused by the common limited spatial extension of the 
radial basis sets, which do not cover the tails of the high energy wavefunctions exceeding 



Rimax- Due to a smaller density of DVR points in these bases, their plateau starts earlier 
than that of V cut = 8.0 eV. Apart from that, the errors of the smaller calculations lie in the 
expected order. Results for the odd states are similar. The fact, that all errors lie below the 
desired precision, proves the convergence of the 3D computation. 

IV. STATISTICAL ANALYSIS OF THE ENERGY LEVELS 



Our result consists of the energies of 1606 even and 1361 odd levels of NO2 |30 
(Note that the symmetries of this adiabatic calculation have no physical meaning for the 
true non-adiabatic molecule with its dominant conical intersection). At energies below 
10 000cm -1 they agree with the published list of experimental values [31J to within an aver- 
age error of 14 cm -1 . Since in that energy region the potential fit has been adjusted using 
the experimental data, this shows the correctness of our calculation in this energy range. 



Similar agreement exists between our computed data and those computed in ref. JT2| . 
At much higher energies a comparison with experimental data is questionable due to the in- 
accuracies of the extrapolated PESs, which we had to use, and the adiabatic approximation. 
However, statistical analyses are meaningful. So, we have performed separate statistical 
analyses of the 1606 even and 1361 odd levels. One expects a transition from regular to 
irregular dynamics with rising energy not later than the conical intersection at 1.289 eV. 
Since only about 75 even and 50 odd states lie below this energy, nearest neighbor histograms 
trying to show this effect, suffer from insufficient data. Nevertheless, in Fig. |2| the nearest 
neighbor distributions of the even states are shown for four different energy windows cov- 
ering all parts of the spectrum. All distributions imply chaotic dynamics by their Wigner 
type shape, but only the two high energy distributions contain enough data to be really 
trustworthy. One may, perhaps, wonder why the adiabatic, uncoupled, ground state PES 
of NO2 shows " chaos" at such a low energy. However, one should not forget that (a) any 
anharmonic PES shows chaos gradually increasing with energy, and (b) that the low-lying 
conical intersection in N0 2 indirectly distorts the shape of the adiabatic surfaces very much 



8 



compared with, e.g., a triple Morse potential. The contribution of the Born-Oppenheimer 
breakdown, i.e. the direct coupling between the surfaces as such, does not change the level 



statistics too much 18 



The rigidity of the spectrum, A 3 (L) seems to be more robust. It is defined as the local mean 
square deviation of the level staircase N(E) from the best fitting straight line over an energy 
range [Ei, Ei + i] corresponding to L mean level spacings, namely 

A 3 (L) =gj5 I j f Ei+L [N(E) -aE- b] 2 dE\ . (10) 

We have calculated it for lengths L between 2 and 300 levels for the whole spectrum, and for 
several fractions of it. We have also used various unfolding procedures |32|] (and no unfolding 
at all, which has no effect at energies above 2.3 eV). We find no energy range in which A 3 (L) 
is near L/15, the expectation value for random level positions, which one generally correlates 
with regular classical dynamics. Instead, taking the usual values of L from 2 to 30, we find 
A 3 (L) slightly above the expectation for a Gaussian orthogonal ensemble (GOE), but never 
higher than twice that value. Interestingly, when we looked for the fluctuations at a much 
wider scale, taking L between 100 and 300, we found an energy dependent behavior, which, 
in addition, differed for the even and odd states. 

An example is shown in Fig. §. For L = 100, the expected values for random and GOE 
sequences are 6.67 and 0.46, respectively. Our results lie between 0.5 and 2.5 in different parts 
of the spectrum. They are robust against a variation of L between 100 and 300, and against 
random shifts of the energy levels of up to 0.5%, simulating errors in the diagonalization 
procedure. To demonstrate that Fig. [| indeed shows the developement of the dynamical 
structure of the states with energy, we show the more conventional way of computing the 
averaged statistics (A 3 (L)) in Fig. [| for two energy ranges: < E < 1.3eV below the 
conical intersection, and 2.5 < E < 2.8eV, where the difference between the symmetries 
is most prominent. The averaged statistics confirms the data of Fig. |3] as representative. 
Therefore, we are convinced that these changes of the rigidity of the spectrum with energy 
are a property of the spectrum of NO2 , or at least of this model of it. 

9 



The saturation of the curves in Fig. |]b was explained by Berry |33| in a semiclassical 



treatment of an integrable and a totally chaotic system. According to this paper, the value 
of L max at which saturation takes place is associated with the period of the shortest classical 
orbit of the system. At this moment, it is unclear to us, why L max is different for even and 
odd symmetries. 

In conclusion, we have shown that the calculation of all J=0-states of a non-hydride second 
row molecule having a few thousand states is feasible, and we foresee that the same holds 
in the presence of a conical intersection, and for low lying resonances. No supercomputer is 
needed if one is willing to pay the price of some patience. 
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FIGURES 
FIG. 1. The reciprocal of the median (taken from a window of 100 eigenvalues) of the density 

of even states p(E), and the median of the absolute errors of the energy eigenvalues of the various 

computations with respect to the largest calculation (with V cu t = 8.1 eV) vs. energy. 

FIG. 2. Nearest neighbor distributions of the even states for four different energy ranges: (a) 
below the conical intersection, (b) about the conical intersection, (c) slightly, and (d) high above 
the conical intersection. The best fits to Poisson and Wigner distributions are also shown. 

FIG. 3. A3 statistics for a window of 100 even (solid line) and odd (dashed line) levels shifted 
in steps of 10 over the energy range. The A3 value is assigned to the energy in the middle of 
the window. The horizontal line is the value of the corresponding A3 for a Gaussian orthogonal 
ensemble (GOE). 

FIG. 4. The averaged A3 statistics and its rms error for a window of L states moved between 
(a) and 1.3 eV, and (b) 2.5 and 2.8 eV. The expectation values for Poisson statistics (P) and the 
Gaussian orthogonal ensemble (GOE) are also shown. 
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TABLES 
TABLE I. Values of the parameters a and b of the incomplete Gamma function used to switch 

between potential parts. 

correction a b [ao 1 ] 

short-range V\\ 1.1 -30.0 

long-range V n 3.5 6.34902 

long-range V 22 2.0 2.09979 

angle V 12 2.0 20.0 
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